Energy flows in vibrated granular media 

Sean McNamara*-^'^^ and Stefan Luding^^^ 
(1) Institute for Computer Applications 1, Pfaffenwaldring 21, 70569 Stuttgart, GERMANY 
(2) Benjamin Levich Institute and Department of Physics, The City College of the City 
University of New York New York, NY 10031, USA 
(February 1, 2008) 

We study vibrated granular media, investigating each of the three com- 
ponents of the energy flow: particle-particle dissipation, energy input at the 
vibrating wall, and particle-wall dissipation. Energy dissipated by interparti- 
cle collisions is well estimated by existing theories when the granular material 
is dilute, and these theories are extended to include rotational kinetic energy. 
When the granular material is dense, the observed particle-particle dissipa- 
tion rate decreases to as little as 2/5 of the theoretical prediction. We observe 
that the rate of energy input is the weight of the granular material times an 
average vibration velocity times a function of the ratio of particle to vibration 
velocity. 'Particle-wall' dissipation has been neglected in all theories up to 
now, but can play an important role when the granular material is dilute. 
The ratio between gravitational potential energy and kinetic energy can vary 
by as much as a factor of 3. Previous simulations and experiments have shown 
that oc y^, with 6 = 2 for dilute granular material, and 5 ~ 1.5 for dense 
granular material. We relate this change in exponent to the departure of 
particle-particle dissipation from its theoretical value. 
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I. INTRODUCTION 



In this paper, we study granular materials where energy is provided by a vibrating plate 
[see Fig. 1(a)]. 




(A) 
(B) 
(C) 



(a) (b) 
FIG. 1. (a) A sketch of the studied system, (b) The wave forms used to drive the vibrating 
plate. 

In this system, two general classes of motion can be imagined, which we will name "coherent" 
and "incoherent" . In the first case, the particles move together in a coherent layer, bouncing 
at some frequency related to the vibration frequency. This state exhibits surface waves, and 
has been the subject of much recent work Another possibility is that the motion 

is incoherent. The particles remain suspended above the vibrating plate, with a density 
profile that - excluding fluctuations - remains constant throughout the vibration period. 
This second state has also been the subject of recent experiments, simulations and theories 

In this paper, we study this second type of motion. We search for a basic understanding 
of the physical processes which govern the input and dissipation of energy. This is a question 
of theoretical interest, because it provides a good test of kinetic theories of granular media. 
These theories, modeled after the kinetic theory of gases ]TT|-pr5[ have been shown to have 



quantitative success only for unforced granular media in the absence of gravity [12,14-16 



However, incoherent vibrated granular materials may be an experimentally accessible system 
well described by these theories. [[7|,p|,p!3| 

The paper is organized as follows. In Sec. |I|, we describe the studied system in detail 
and set forth our notation. In Sec. |IT1|, we review the literature. In Sec. |rV|, we check the 
previous results against our simulations, and construct some new theories. In particular, 
we examine five topics: 1) the rate of energy dissipation in particle-particle collisions, 2) 
the effect of particle rotations, 3) the energy input by the vibrating floor, 4) the energy 
dissipated by particle-wall collisions, and 5) the ratio of kinetic to gravitational potential 
energy. FinaUy, in Sec. 0, we assemble the best theories and investigate the dependence of 
the energy on vibration velocity. 
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II. DESCRIPTION OF THE SYSTEM 



A. Definitions 

A sketch of the system is shown in Fig. |l](a). A granular medium, modeled by a gas of 
inelastic disks, is contained in a box of width L and infinite height. The gas is pushed against 
the bottom by a gravitational field with acceleration g. Energy is added to the system by 
the bottom of the box, which vibrates with period r and typical velocity V . (In Sec. [IV (J| , 
we relate V to the maximum velocity Knax-) The granular medium consists of particles 
of radius a and mass m. Convenient nondimensional parameters describing the system are 
the length measured in particle radii, the number of layers of particles at rest H = 2aN/L, 
and the number of particle radii a particle, initially at rest, falls during one cycle of a wall 
vibration, gT'^/{2a). The description of the system is completed by specifying the particle- 
particle and the particle-wall collisions, and the wave form of the vibrating wall. We use the 
three wave forms shown in Fig. ^(b). Wave forms A and B are convenient for theoretical 
analysis and wave form C is a computational convenient numerical approximation 



to a sine wave, constructed by patching together parabolas. In Sec. [IV (J| , we will estimate 
the effect of this approximation on the energy input at the vibrating wall. 

We consider collisions with constant normal and tangential restitution coefficients. This 
collision model has been widely studied and is relatively easy to analyze theoretically. The 
normal (tangential) restitution coefficient for collisions between particles is denoted Vp {(3p), 
and for coUisions between the particles and the side walls (/?«,)• The vibrating floor is 
elastic (r = 1) and smooth {(3 = —1). More details on the collision rule are presented in 
Appendix ^. 

We describe the state of the system with three types of energy, each defined per degree of 

_ o 

freedom, per particle: E, the translational kinetic energy of the particles, E, their rotational 
energy, and Eg, their gravitational potential energy. These energies are calculated as follows: 

N 2 N -\ N 

- m ^ ° mqa ^ n ^ 1 x-^/ , x 

^-Im^--' ^ = = "9]^ gfe (1) 

Here, Vi is the translational velocity of the i^^ particle, Ui is its angular velocity, and Ui is 
its height above the time averaged position of the vibrating floor. The center of mass of the 
particles at rest is ho. We calculate ho from Eq. (14) of Ref. |^. The number of translational 
(rotational) degrees of freedom per particle is n {n). In this paper, we consider exclusively 

_ o 

two dimensional systems, where n = 2 and n = 1. E and E are often called "granular 
temperatures" . This terminology is not meant to imply that a thermal equilibrium exists in 
granular flows, but simply to draw an analogy between these quantities and the temperature 
of an ideal gas, which is proportional to the average energy per degree of freedom. 



B. General Approach 

_ o 

The "holy grail" of this paper are expressions for the energies E, E, and Eg in terms of 
the parameters L, g, V, a, r, m, Vp, r^, Pp and together with a physical understanding 
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of the system. The starting point of the quest is the equation for the energy flow at steady 
state: 



Pb — Dpp ~l" Dpwi (2) 

where Ph is the power input by the bottom, Dpp is the power dissipated by particle-particle 
collisions, and Dp^ is the power dissipated by particle-wall collisions. We will try to express 
each of these quantities in terms of E and the system parameters. Then Eq. (Q) will give 
an expression for and we will also understand what processes determine E. 

In simulations, it is possible to monitor each of the three quantities in Eq. @ or set each 
of them either to zero or to a known constant. Dpp can be set to by using elastic particles, 
just as using elastic walls permits one to set Dp^ = 0. Pb can be set to a known value by 
driving the bottom wall with waveform B [0. We will discuss each of the three terms in 



Eq. (ID and exploit this separability to construct and verify our theories piece by piece. 



C. Simulational Approach 

We use the standard event-driven method, which has been previously used to study 
vibrated granular materials [p|,^ p!8| , p!9| . Its main approximation is to consider that collisions 



are instantaneous, i.e. the time of contact during a collision is 0. The main disadvantage of 
event-driven simulations for granular materials is the occurrence of inelastic collapse - an 
infinite number of collisions occurring in a finite time pO| , pT| . To circumvent this problem, 
dissipation is turned off when the inelastic collapse singularity is approached. A similar 
method has been used in other recent work P,p^. The fraction of dissipationless collisions 
is always less than 10~^ for the simulations shown here. In addition, we repeated the 
three simulations that had the most dissipationless collisions, approaching more closely the 
inelastic collapse singularity, and found no measurable change 

In the simulations presented here, we use the particle radius to define the unit of distance, 
and the particle mass to define the unit of mass. The unit of time is arbitrary. To apply 
our results to experimental systems, it is necessary to introduce conversion factors based on 
the particle mass and the particle radius, and to choose r and g so that gr^/a has the same 
value as the experimental system. 

In order to vary the system parameters in a consistent and organized way, the following 
procedure was used. First, 'central' values were selected for all the parameters except V. 
Then V is swept from small values to very large values, care being taken that all simulations 
fall within the "incoherent" category discussed above. The central values used in this paper 
are = 160, L = 50, Vp = 0.95 = 1, g = 1.0 and r = 1.0. Then for each parameter, 
two series of simulations are run, one where the parameter is increased by a factor of 5, and 
another where it is decreased by a factor of 5. There are two exceptions: 1 — Vp (not Vp) is 
increased or decreased by a factor of 5, and g is increased or decreased by a factor of 25. 

Most of the quantities measured during the simulations must be averaged over long 
periods of time to obtain stable results. To obtain a suitable averaging time, an "energy 
turn-over time" te = E/Pb was estimated. The averaging time was taken to be 20te- It 
was verified that the actual te was always close to the estimated one. Initial transients by 
running each simulation over several averaging periods, and verifying that averages were not 
changing significantly with time. 
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III. REVIEW OF PREVIOUS WORK 



A. Experiments and theory of Warr, et al. 

Warr, Huntley and Jacques modeled vibrated granular media as an isothermal "at- 
mosphere" . By integrating over the "atmosphere" , one obtains 

= 2(1 - rl)NmgH{E/mYl\ (3) 

Next, one can estimate the energy added by considering the density at the bottom of the 
"atmosphere" . The result for wave form A is 

Pfc = {l/2)NmgV'^{m/Ef/'^, (4) 

Setting Pf, = Dpp yields 

mV'^ . . 

^ = 4i7(l-r2)- 

This result was calculated for smooth particles = — 1), with the system being driven by 
wave form A. The authors also present experimental results suggesting 

E^E-^'V'^ Eg mgH~''^V^<> . (6) 

with z/ = 0.6 ± 0.03, 6 = 1.41 ± 0.03, Ug = 0.27 ± 0.11 and 5^ = 1.3 ± 0.04. In contrast, the 
theory assumes that the gravitational potential energy is always proportional to the kinetic 
energy with u = Ug = 1 and 6 = Sg = 2. Explaining the discrepancy between Eq. (§) and 
Eq. (0) was one of the major motivations for our study. 



B. Simulations of Luding 

Luding and co-workers have studied this problem with numerical simulations. Simula- 
tions without rotation [|1^] give Ug ^ 1, 6g ~ 1.5. Simulations including rotation, carefully 
planned to duplicate the experiments give Ug = 0.76 ±0.11 and Sg = 1.60 ±0.10. When V 
is made very large and dissipation at the walls is suppressed P,pi]|, 6^2. Lee |^ obtained 



results that are qualitatively similar. Luding [|1^] used a shghtly different collision model 
which accounts for Coulomb friction. In the limit of infinitely strong friction, it reduces to 
the model used in this paper. 



C. Theory of Kumaran 

Kumaran followed the same approach as Warr, Huntley and Jacques [0, except that 
he assumes that the particle velocities are distributed according to a Maxwellian velocity 
distribution. On the other hand, Warr, et. al assumed that all particles possess the mean 
velocity. Thus, Kumaran obtains the same results as Warr, et. al, except for the constant 
prefactor. Kumaran also investigates both wave forms A and B. He obtains 
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Dpp = J'R/2{l-rl)NmgH{E/mf'^, and (7a) 



n = Nmg [{V) + 2y^(V2) + 0{mV^)^ , (7b) 

where the angle brackets indicate an average taken over one period. If a wave form is 
symmetric, hke wave forms A and C, then {V) = 0, and Pb ~ NmgV'^{m/Ey/'^, as in 
Eq. d^), but when the wave form is asymmetric (wave form B), then the first term of the 
series is nonzero and dominant, so that Pf, ~ NmgV. 



D. Theory and Simulations of McNamara and Barrat 



McNamara and Barrat JT^ studied the input of energy at a vibrating wall in the absence 



of gravity. They also found a difference between wave form B and wave forms A and C. For 
wave form B, the power input by the vibration wall is Pi, = pVL, where p is the average 
pressure on the vibrating wall. This equation can be derived by considering the encounter 
of a single particle with a moving wall at time t^^. The particle's change in energy AE is 
related to its change in momentum Ap by AE = V{t^:)Ap, where V{t^) is the wall velocity 
at time t*. For wave form B, V{t^,) = V always, so it is easy to average over time and find 
that the energy input is pVL. For wave forms A and C, the velocity of the plate is not 
always the same when the particles hit the plate. However, the probability of a particle 
hitting the plate at a given phase is governed only by the quantity U/V, where U is the 
velocity of the incoming particle. Therefore, for wave forms A and C, Pb = pVLf{U/V), 
where the function / is unknown. To adopt these results to a system under gravity, one sets 
pL = Nmg, the weight of the granular material. The result for wave form B, 

Pb = NmgV, (8) 

agrees with Eq. (|7^) for asymmetric wave forms. For wave forms A and C, Pb = 
NmgV f (U/V), which coincides with Eqs. (|) and (|7bD if fiU/V) = V/U. 



IV. RESULTS 
A. Particle- Particle dissipation of smooth particles 

1. Simulations 

The two expressions for Dpp in Eq. (|^) and Eq. ( |7aD differ only by a constant. We measure 
this constant by calculating 

f = ^PP (n) 

{l~rp)NmgH{E/my/^' ^ ' 

Eq. (^ predicts Cpp = 4 and Eq. ( |7aD predicts Cpp = \phx ~ 2.5. [The factor of 2 difference 
between the predicted values Cpp and Eqs. (0) and ([7a| ) arises from replacing \ — ri with 



6 



2(1 — Tp). For Tp ^ 1, 1 — Tp ^ 2(1 — Tp), so one might think that these two quantities would 
be equivalent. However, the difference at = 0.75 is sufficiently great to show that 1 — 
collapses the data better than 1 — r^.] Fig. ^ shows that the scaling successfully collapses 
the simulation data onto a single curve. 



O 




10"' 10° 10' 10' 10^ 10^ 

(h-h„)/a 

FIG. 2. A test of the scalings in Eqs. (^) and (|7a|). Here, Cpp is plotted against the height of the 
center of mass h = J2 Vi above its position at zero energy, /iq. Both equations predict that Cpp 
should be a constant, which is indeed true for large heights (dilute granular media). The central 
simulation has N = 160, rp = 0.95, L = 50, 5 = 1, and r = 1; ^ is varied over several orders of 
magnitude. The other series of simulations, each indicated by a different symbol, have the same 
parameters as the central one, except for the parameter shown in the legend. The bottom is driven 
with wave form B, and the side walls are replaced by periodic boundary conditions. The dashed 
line shows the constant predicted by Eq. ([7a]). 

The figure can be divided into two regions: the dilute region, h — > 50a, and the dense 
region, h — < 50a. In the dilute region, Cpp is a constant, and quite close to the value 
predicted by Eq. (^). The points collapse into families determined by the particular value 
of if = 2aN/L of each simulation. However, the scaling captures the dependence on V, 
Tp, g and r to within the noise in the simulations. In the dense region, h — < 50a, 
all simulations collapse tightly onto a curve Cpp ~ \og[{h — ho) /a]. The reason for this 
dependence of Cpp on {h — hQ)/a is unknown. It is not due to the increasing density, because 
the low H {N = 32 and L = 250 curves in Fig. 0), as well as additional simulations (not 
shown) also collapse onto the same curve. Nor is the deviation due to waves propagating 
upwards from the vibrating plate; the agreemt of the r = 0.2 curve with the rest of the 
simulations excludes this. We also tried an alternate way of adding energy: the bottom 
plate is held fixed, and when a particle hits the bottom, it is given a velocity drawn from 
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a Maxwellian distribution. But these simulations also reproduce the curve shown in Fig. |^. 
Replacing periodic boundaries with elastic walls modifies only slightly Fig. ^ 

Fitting a straight line to the points with a < h — ho < 50a (and excluding the = 800 
and L = 10 points) gives Cpp ~ 0.301og[(/i — /io)/o] + 1-35. This relation is valid over 
only about one and half orders of magnitude, so it is difficult to say whether it is really 
logarithmic, or whether the logarithm is only a convenient approximation. 



2. A physical argument 

We now derive the scaling relation from physical arguments. This permits us not only 
to understand why 1 — Tp is better than l — r^, but also enables us to extend the theories to 
account for particle rotations and dissipation at the side walls. 



From general kinetic theory arguments in the style of Haff , we argue that the dissi- 
pation due to collisions between particles will be 



i^E) Q , (10) 



where (AE) is the average energy lost per collision, f/ is a typical particle velocity, and s 
is a typical particle separation. Multiplying (AE) by the collision frequency U/s gives the 
energy dissipation per particle per unit time. Then multiplying by the number of particles 
gives the total energy dissipation in the system. 

Next, we must relate (AE), U, and s to ^ and the independent parameters. First of 
all, Eq. (|A9|) in the appendix shows that (AE) ~ (1 — rp)E for smooth (/3p = —1) particles. 
Next, since U is a typical velocity, we make the identification U ~ {E/niY^'^. Finally, to 
estimate the particle separation s, we turn to a one dimensional model. We imagine a column 
of n particles, suspended in a gravitational field by the vibrating plate at the bottom of the 
column. The column is in a steady state, so the net force due to collisions on each particle, 
averaged over time, must balance the gravitational force mg. Consider a particle somewhere 
in the column, with M particles above it (1 < M < ra). In a steady state, this particle 
must receive a momentum fiux (M + l)mg from the particle below it and transfer Mmg 
to the particle above it. The momentum transferred per collision is (1 + rp)mU, and the 
collision frequency is U/s, hence the momentum fiux will be proportional to m(l +rp)U'^/s. 
Equating the two expressions for the momentum fiux 

Mmg ~ (1 + rp)mU^/s, (11) 

and solving for s gives s ~ (1 + rp)f/^/ (Mg). M will scale as n, and in two dimensions n can 
be approximated by the number of layers of particles H = 2Na/L. Again using t/^ ~ E/m, 
we have s ~ (1 + rp)E /{mgH). 

Inserting our expressions for s and U into Eq. (|ToD , we have 

Dpp = Cpp{l - rp)NmgH{E/mf/\ (12) 

This agrees in order of magnitude with the scaling tested in Fig. |^, and, in addition, explains 
why the 1 — rp of Eq. (|^) collapses the data better than the 1 — of Eqs. (0) and (|7aD. 
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B. Particle-particle dissipation for rough particles 



All the simulations and theories presented so far in this paper consider "smooth" particles 
{Pp = —1), thus ignoring rotation. But particles in granular flows rotate, and our theory 
would not be complete without considering rotation. 

We now revise our argument for the scaling of Dpp in Sec. [IV A| , to take into account 
rotation. The only place the smoothness of the particles entered was in estimating the 
change in energy per collision (AE) in Eq. ([T0|), so we use the more complicated expression 
given in the appendix, Eq. (|A9|) : 



(AE) = -2(1 - rl)E 



g(i-/3p^) 

l + q 



E + E/q 



(13) 



This involves E as well as E, so we draw from a recent paper on the ratio of translational 
to rotational kinetic energy in granular flows. One of the principle results of this paper is 
that 



E 

K=- 
E 



l + 2q-(3p 

q{i + Pp) 



(14) 



in vibrated granular flows. This convenient fact enables us to replace E with E/K. The 
result is 



{AE} 



E 



-2[1 -r;']E. 



l + 2q- j3p_ 

with the "effective restitution coefficient" , governing the loss of energy. 



.2 - 



l + 2q-(3p 



Inserting {AE) 



)E into Eq. (|T0|) gives 



'1 -r! 



n — r'° 



NmgH{E/m 



,1/2 



(15) 



(16) 



(17) 



as the analogy of Eq. (|1^). This equation also defines a constant C°p, in analogy with Cj 
We plot C;^ in Fig. |. 



pp- 
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(h-h„)/a 

FIG. 3. A test of the scaling in Eq. ([T7|) , which defines C°p. The parameters for the simulations 
are N = 160, L = 50, 5 = 1, r = 1, and Vp = 0.95, with (3p given in the figure and V varied between 
100 and 0.35. The dotted line and the /3p = — 1 simulation also appear in Fig. |2[ Wave form B is 
used, and the side walls are replaced by periodic boundaries. 

There remains a weak dependence on /3p, but Cpp = Cpp within the accuracy of the simula- 
tions. 



C. Energy input 

We now turn our attention to the energy added by the vibrating floor. We test first 
the simplest result: Eq. (^), which applies only to wave form B. Reviewing all simulations 
shown in Fig. |^, we find the largest deviation from Eq. (H) is 1.6%, with most others much 
less (90% of the particles deviate by 0.5% or less). This equation is accurate because it 
is independent of the velocity distribution of the particles, and can be derived from the 
conservation of momentum |17| ]. 

For wave form A, Refs. |^J^ predict ~ NmgV^{m/Ey^'^. We test this scaling against 
the simulations in Fig. ^. 
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'10"' 10° 10' 10' 10^ 10^ 

(h-h„)/a 

FIG. 4. A test of Eqs. (|) and @: we plot Cb = Pb/[NmgV'^{m/ E)^/'^] for the symmetric 
wave form A. The symbols and parameters of the simulations are the same as in Fig. ^, except 
the bottom plate is driven by wave form A instead of wave form B. The dashed lines show the 
constants predicted by Eqs. (§) and ([7b|). 

We see that the scahng is only partially successful. The rescaled power input varies by as 
much as a factor of 5, and there is a strong dependence on H{1 — Vp) not captured by the 
scaling. An analogous plot for wave form C is similar. 

To find a better way to calculate for wave forms A and C, we try the scaling suggested 
by Ref. = NmgVf{U/V) with U = {E/my/^ Let us now consider the unknown 

function f{U /V) in the limits U/V — > oo and U/V —>■ 0. In the first limit, the particles are 
moving infinitely more quickly than the wall, and thus have an equal probability of hitting 
the wall when it is ascending or descending. In this case, the net energy input by the wall 
will be because the energy gained by particles during the ascending phase is lost during the 
descending phase. As U becomes smaller, the particles have a higher probability to hit the 
wall during its ascending phase, and hence Pb becomes positive. Finally, in the limit ?7 — > 0, 
the particles can hit the wall only during the ascending phase, so wave form A becomes 
equivalent to wave form B, so f{U/V) 1. This limiting value of / is also attained for 
wave form C if V = 2Vmax/3. In all cases, V is equal to a space average of the velocity of 
the plate during its ascending phase, V = A~^ J^V{y)dy, where y = A is the maximum 
height attained by the plate. When the plate sweeps through a motionless gas of particles, 
V is the average plate velocity seen by the particles (assuming all particles collide only once 
with the plate). For the sinusoidal wave form, a similar calculation gives V = nVraa.x/4:. 

The unknown function / is shown in Fig. |^ for wave form A. An analogous plot for wave 
form C is similar. 
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FIG. 5. A graph of f{U/V) = Pb/ {NmgV), with U = {E/m)^^'^. The parameters and symbols 
of the simulations are the same as in Fig. |2|, except that wave form A is used. 

For wave forms A and C, points fall on a single curve, but they are bunched together 
depending on the value of {l—rp)H. The function / is well approximated by an exponential, 
so that 

Ph = NmgVexp[-aU/V]. (18) 

Least squares fits give = 0.353 ± 0.002, and ac = 0.475 ± 0.003. (Of course, = 0.) 
We note that a^i/ac ~ 3/4. Fig. ^does not prove that / is an exponential, because we have 
calculated / over less than one order of magnitude. However, an exponential is a better 
approximation for / than a power law or a straight line. 

Note that in all scalings, P}, depends on the plate motion only through V and not through 
r. This differs from previous simulations, which suggest that Pi, drops dramatically when r 
is below a critical value [Q. However, Ref. @] used a soft sphere simulation method, where 
particles remain in contact during a finite time while colliding. If r approaches this contact 
time, the efficiency of the forcing will drop. On the other hand, if r is made very large, the 
granular material dissipates most of its energy during one vibration cycle, and the transition 
to the "coherent" state occurs. In the coherent state, the relevant parameter describing the 
forcing is no longer V, but the acceleration V/t. 



D. Effect of Side Walls 

Finally, we consider dissipation of energy through collisions between the particles and 
the side walls. There is no theory for Dp^j in the literature. Nevertheless, we find that Dp^ 
can be estimated by 
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= CpMDpp (I) l^^j , (19) 

where v = {h — ho)/{aH) = L(h — hQ)/{2Na?) is proportional to the specific volume or 
inverse density. Cpwiv) is a linear function, and is the kinetic energy per degree of 
freedom in the x component of the particle velocities: Ex = 'm{2N)~^ J^v^i- In a gas at 
thermal equilibrium, E^/E = 1, but in these simulations, E^/E can attain values as high 
as 0.2. The function Cpwiv) is shown with the data from the simulations in Fig. ^. 
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FIG. 6. A test of the scaling in Eq. ([T9|), where Cpw = [Dp^ / Dpp){L / a){E / E^f'/'^ is linear in 
V = {h-ho)/{aH) = L{h-ho)/{2Na'^). The dotted line is Cp^ = 1.45 + 0. 905t;. These values were 
obtained by fitting all the points on the graph except those of the large H simulations (L = 10 and 

= 800). These simulations have the same parameter values as those in Figs. ^, ^ and|5|, except 



that periodic boundaries have been replaced by walls with r«, = Vp and /3„ 



-1. 



The simulations agree well with Eq. ([T9|) , except for the large H (L = 10 and = 800) 
simulations. 

The result Eq. (^) can be understood as a continuation between two limits. In the limit 
of t; < 1, Eq. (ig) becomes Dp^ ~ Dpp{a/L){Ex/Ef/\ The factor of {E^/Ef/^ arise the 
X velocities alone determine the frequency of collisions with the wall, and the amount of 
energy lost. The factor oi a/L appears because in this limit, the particles remain tightly 
packed, and rarely change places. The ratio between Dpp and Dp^ will be equal to the ratio 
between the number of particle-particle contacts and the number of particle-wall contacts. 
Geometrical considerations show that these ratios are of order a/L. 

In the limit oi v ^ 1, Eq. ([T^) becomes [after using Eq.(|7aD for Dp^] Dp^ ~ (1 — 
rp)A^£'|/^/(m^/^L). Then, realizing that the average energy lost in particle wall collisions 
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is {Ew) ~ (1 — r1j)Ex (note that = rp), and that the typical x velocity will be Ux = 
{Ex/mf^, Eq. (|l|) becomes 

Dp^r^N{E^){UjL). (20) 

The similarity between this equation and Eq. ( p!0D for Dpp permits the following interpreta- 
tion: in the large h limit, the particles are bouncing back and forth between the two side 
walls. The factor U^/L gives the frequency of collisions with the wall, so that {Ew){Uxl L) 
gives the rate of energy loss per particle. There are N particles, so a factor of appears in 
Eq. (li. 

The deviation of the large H simulations from Eq. ( [19|) is probably due to their different 
mass distribution. These simulations differ from the others because they form a dense plug of 
particles suspended by a dilute "hot" region. In contrast, the others have a density maximum 
near the bottom, similar to a normal atmosphere. An additional curiousity is convection 
which appears in the large N simulations (solid circles in Fig. ^). Fig. |^ shows the motion 
of all the particles during one cycle for a N = 800 simulation, revealing a circulation of 
particles within the plug. When V is decreased below a critical value, the circulation ceases, 
and this critical value of V corresponds to the break in the N = 800 curve (solid circles) in 
Fig. ^. We have not studied this convection in detail. 




FIG. 7. Streaklines showing the motion of each particle during one cycle of the wall vibration. 
The solid dot shows the position of the particle at the beginning of the wall vibration, and the line 
shows its motion during the cycle. The simulation has = 800, = r^, = 0.95, /?p = = — 1, 
L = 50, g = 1, T = 1, and V = 15.5 with wave form B. The figure is horizontal, the vibrating 
bottom is on the left. 

A further complication arises when the walls are rough [f]^ ^ ~1)- One must decide 
whether the walls move with the vibrating bottom or not; i.e. whether the bottom is a piston 
moving vertically between the stationary walls, or whether the particles are contained in a 
box which is shaken. If the side walls also move, then energy can be input at the side walls 
as well as at the bottom (i.e. Dp^j can be negative) 0. In this paper, we always consider 
the side walls to be stationary. 



E. Gravitational Potential energy 



As shown in Sec. |1V A| , Dpp is a known function of the system parameters, the energy 
and h — ho- We would like to eliminate this dependence on h — ho before combining our 
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results for Dpp, Pb, and Dp^j. In order to do this, we plot in Fig. || the ratio E /{rpEg) [where 
Eg = mg{h — /iq) is the gravitational potential energy per particle]. 




FIG. 8. The ratio of translational kinetic energy E to gravitational kinetic energy Eg scaled by 
1/rp and plotted against v = {h — h(j)/{aH) = L{h — h(j)/{2Na'^), for the simulations in Fig. ^. 

This scaling successfully collapses the data onto a single curve reminiscent of Fig. ^ We note 
that h must be carefully calculated; if h is taken to be the height above the lowest position 
of the floor instead of the height above the average position of the floor, the simulations 
with qt"^ = 25 fall off the curve. 

Gravitational energy and kinetic energy are often assumed to be equivalent. Fig. |^ 
should be a cautionary note. The dependence of E/Eg on the system parameters is quite 
complicated, and has not been theoretically investigated. 



V. SUMMARY AND TEST 

To summarize and test the formulas for Dpp, Pb, and Dp^ presented in the previous 
sections, we calculate the predicted value of E for a set of simulations from a previous paper 
0. Using Eqs. (^7^, (|l5) and (|19D to rewrite the energy balance, Eq. (||), as an equation in 
terms of the system parameters and U = (E/mY^'^, gives 

gVeM-c^U/V) = Cpp^-^gHU + Cpp^-^gHUCp^{v) (^-j i-^j , (21) 

where Cpp = min(0.30 log[/i — ho] + 1.35, V^), and Cp^ = 1.45 + 0.905i;. In writing down 
Eq. (pTD, we have assumed that Eq. ([T6|) holds for the particle- wall restitution coefficients. 



15 



To close Eq. (pT]), expressions for {h — ho) /a and need to be supplied. We assume 
{h — ho) /a = 2E/{mga) and = E, which are true in a gas of elastic hard spheres under 
gravity (r^ = 1, (3p = ±1, V = 0). These assumptions hold approximately for granular 
media, except in some extreme cases. 

We test this equation against simulations in Fig. ^ 




FIG. 9. A test of the theoretical results from Eq. (|2l|). The central simulations have = 50, 
a = 0.05 cm, L = 20a, Vp = 0.9, g = 981 cm/s^, r = 0.01 s, = 1.0, and fip = (3w = —1, with 
wave form C. The other simulations are the same, except for the parameters given in the figure 
legend. The solid lines give the results of the theory presented in Eq. (]2l|). The dotted lines give 
the theory without taking into account the logarithmic dependence of Cpp on h — ho, i.e., with 
Cpp = \/27r. 

The parameters of the central simulations are = 50, a = 0.05 cm, L = 20a, Vp = 0.9, g = 
981 cm/s^, T = 0.01 s, = 1.0, and Pp = Pw = —1, with wave form C. These parameters 
were chosen to mimic the simulations in Fig. 11 of Ref. 0. There remain three differences 
between these simulations and those of Ref. 0. First, Ref. P] uses a more complicated 
collision rule, where the tangential restitution coefficient varies between —1 and an upper 
limit /3o, depending on the impact parameter. In our simulations, the tangential restitution 
coefficient is the same for ah collisions. Secondly, Ref. used slightly polydisperse spheres, 
whereas we use monodisperse spheres. Finally, we replace the sinusoidal wave form of Ref. ^ 
by wave form C. In spite of these differences, the simulations presented here show the same 
behavior as those of Ref. . 

We show two versions of the theory, one which takes into account the logarithmic de- 
pendence of Cpp (the solid hues) and the other which does not (the dotted lines). The solid 
lines reproduce the observed dependence E ~ V^^'^ for \4iax < 100 cm/s, whereas the dotted 
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lines show E ^ V"^. Thus, the puzzhng scahng observed in previous work is due to the 
logarithmic dependence of Cpp on h — Hq. This dependence does not yet have a theoretical 
explanation. 

For Vmax > 100 cm/s, there is a significant disagreement for the simulations with dissi- 
pative walls (r^ 7^ 1), which can be traced back to the failure of the assumption = E. 

VI. CONCLUSIONS 

This paper has studied the two dimensional vibrated granular media, using the energy 
balance Eq. to organize our investigation. Reviewing existing theories, we find the 
particle-particle dissipation is well understood in the dilute limit. We were able to show 
that Eq. ( [7a]) is very accurate in the absence of side walls. This result was extended to deal 
with rough rotating particles with a constant tangential restitution coefficient. More work 
is needed to account for the more realistic situation of variable tangential restitution. In the 
dilute limit, particle-particle dissipation is well accounted for by existing theories, but in the 
dense limit, it shows an unexplained dependence on the height of the center of mass above 
the vibrating bottom, h — ho, well approximated by a logarithm. It is this departure from the 
dilute theory which accounts for the E ~ V^^'^ scaling observed in previous work. The cause 
of this departure is not yet known: it is not due to a change in density, since it depends only 
on the height of the center of mass. Nor can it be explained by waves propagating upwards 
from the bottom. Furthermore, the ratio between potential and kinetic energy, E/Eg, can 
vary by almost a factor of 3. This variation is not understood, and hampers the ability of 
our theory to predict h — ho. 

The energy input by the vibrating plate is also well known for the special driving wave 
form B in Fig. ^(b). For more conventional wave forms, the energy input is well approxi- 
mated by the exponentially decaying function of U /V shown in Eq. ([T8|). 

To our knowledge, this is the first paper to treat particle-wall dissipation in detail. Our 
theory accurately predicts energy losses at the wall, but requires knowing the difference 
between the horizontal and vertical kinetic energies. A more complete theory would have to 
estimate this difference from the system parameters. 

We acknowledge the generous support of the Alexander von Humboldt-Stiftung and the 
DFG, SFB 382 (A6). 

APPENDIX A: THE COLLISION RULE 

Consider a collision between two particles of radius a. Let Vi be the translational velocity 
of the first particle, uji its angular velocity and ri its position at the time of contact. The 
quantities V2, UJ2, and r2 are the analogous quantities for the second particle. Then the 
relative velocity at the point of contact is 

Vc = Vi — V2 — a{uji + UJ2) X n. (Al) 

Here, the unit vector n = (ri — r'2) / |ri — r2\ points along the line connecting the centers 
of the colliding particles, from particle 2 towards particle 1. The change in the normal 
component of is parameterized by the "coefficient of normal restitution" r^, so that 
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v'^ h = ~rp{vc ■ n), where the prime denotes the velocity after the colhsion. When = 1, 
energy is conserved, and energy dissipation requires < < 1. The coefficient of tangential 
restitution Pp is defined analogously, i.e. v[ x h = —[3p{yc x n). Energy is conserved for 
Pp = —1 (perfectly smooth surfaces) and for Pp = 1 (perfectly rough surfaces). In the ffist 
case the a;, have no effect on outcome of the collision, and do not change during the collision. 
Energy is dissipated when Pp lies between these two extremes. 

From the definitions of Vp and Pp and the assumption that the interaction takes place 
only at the point of contact, it is possible to derive the collision rules 



, l + Tp q{l + Ppj 

^1 2 = ^1,2 T ^ Vn T ^ , ^ [Vt + Vr), and 
2 2q + 2 

where the equation for particle 1 (2) takes the — (+) sign in the top line above. To derive 
Eq. (^^) we used momentum conservation and the definitions 

Vn = [(■^i - V2) ■ n] ■ n, 
Vt = Vi - V2 - Vn, and 

Vr = —a(u;i + 0^2) X n. (A3) 

Here, Vn is the normal component of Vc, Vt is the tangential component due to translation, 
and Vr is the tangential component due to rotation. Note that v^ = Vn + Vt + Vr- 
The change in translational energy is 

= -Qvl - S [Ctiv^ + Ct2{vt ■ Vr) - Ctsv^] , (A4) 

with the positive prefactors Q = m(l — r^)/4, S = mg(l + /3p)/[4(l + g)^], and the constants 
Cti = 2 + g(l — Pp), Ct2 = 2 — 2qPp and Cts = g(l + Pp). Likewise, the change in rotational 
energy is 



AE = -S 



CrlV^ + Cr2{Vt ■ Vr) + CrsV^r . (A5) 



where the constants are Cri = (1 + Pp), Cr2 = 2(g — Pp), and Cr3 = 2q + 1 — Pp. Note that 
the C are positive (only Cr2 can also be negative) so that the signs in Eqs. (|A4| ) and ( |A5|) 
indicate the direction of energy transfer between the degrees of freedom. 
Eqs. ( |A4|) and (|A5|) can be added together to give 

AE = -Qvl - S{1 + g)(l - Pp){vt + Vr)\ (A6) 

In this paper, we need to know the average energy lost per collision. Using angle brackets 
to denote averages over collisions, we have 

{AE) = -Q{vl) - S{1 + g)(l - Pp){{vt + Vr)'). (A7) 

Assuming that the particles' velocities are distributed according to a Maxwellian velocity 
distribution, and that their postions and velocities are uncorrelated gives 

{vD = 8E/m, {{vt + VrY) = 4(E + E)/m, (A8) 
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in two dimensions. ||23| Thus the total energy lost during one collision is 

9(1 - Pi) 



{AE) 



-2{1 -r;)E 



E + E/q 



(A9) 
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